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SUMMARY 


Power  series  expansions  (with  coefficients  obtained  by  recurrence 
formulas)  are  more  efficient  than  other  integration  procedures  for 
computing  concurrently  an  orbit  and  the  resolvent  matrix  of  its  varia¬ 
tional  equations,  in  the  Restricted  Problem  of  Three  Bodies.  For  the 
same  requirements  on  accuracy,  the  series  expansions  use  only  about 
30  per  cent  of  the  computing  time  of  the  multi-step  procedures,  and 
only  12  to  15  per  cent  of  the  computing  time  of  the  Runge-Kutta-Nystrom 
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I .  INTRODUCTION 

Steffensen  [1]^  has  set  up  recurrence  formulas  to  integrate  by 
power  series  the  planar  Restricted  Problem  of  Three  Bodies;  the  equations 
of  motion  were  taken  in  their  Lagrangian  form,  in  the  jovicentric 
synodical  coordinate  system.  The  series  involved  were  proved  to  be 
convergent  as  long  as  the  initial  point  is  not  lying  at  one  of  the 
singularities.  Steffensen' s  algorithm  has  been  used  for  the  first  time, 
and  quite  extensively,  by  Rabe  in  determining  the  long  period  orbits 
around  the  Lagrangian  equilateral  centers  of  libration,  either  in  the 
Sun-Jupiter  system  [2]  [3]  or  in  the  Earth-Moon  case  [4J.  Fehlberg  [5] 
adapted  the  method  to  the  Lagrangian  equations  of  the  Restricted  Problem 
in  the  barycentric  synodical  coordinate  system;  he  also  extended  it  to 
the  problem  of  a  charged  particle  in  the  field  of  a  magnetic  dipole.  In 
both  problems,  integration  by  recurrent  power  series  proved  to  be  faster 
and  more  reliable  than  other  integration  procedures. 

Recently  Deprit  and  Price  [6]  have  brought  Steffensen's  ideas  into 
play  to  integrate  concurrently  the  Hamiltonian  equations  of  the  Restricted 
Problem  and  their  related  variational  equations;  in  particular,  they 
propose  to  integrate  in  one  package  a  periodic  orbit,  the  4  *  4-matrix 
of  its  fundamental  displacements,  and  its  characteristic  exponents. 

It  amounts  to  repeatedly  evaluating  by  recurrence  33  power  series;  the 
time  step  is  variable,  being  adjusted  automatically  by  the  program  anywhere 
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1)  Numbers  in  square  brackets  refer  to  References,  page  14. 
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within  a  preassigned  interval.  Error  control  is  exerted  by  5  first 
integrals  and  6  bilinear  identities  derived  from  the  symplectic  character 
of  the  resolvent  matrix. 

Numerical  integration  by  recurrent  power  series  is  not  a  general 
purpose  algorithm:  in  each  case,  the  right  hand  members  of  the  differen¬ 
tial  system  should  be  manipulated,  and  the  order  itself  of  the  system 
possibly  increased,  so  as  to  give  to  the  equations  the  general  form 
considered  by  Fehlberg.  But  such  an  ad  hoc  treatment  has  precisely  the 
virtue  of  suiting  the  numerical  integration  method  to  the  particular 
problem  under  investigation.  It  often  means  greater  reliability  on  very 
long  spans  of  integration,  unusually  large  step-sizes,  and  an  appreciable 
saving  in  computing  time. 

We  checked  these  benefits  of  the  power  series  in  the  complete 
(orbit  +  variations)  integration  of  the  Restricted  Problem  against  more 
classical  methods,  like  the  Runge-Kutta-Nystrom  algorithm  or  the  multi- 
step  procedures.  To  reach  an  accuracy  of  at  least  9  decimal  figures  for 
the  characteristic  exponents  of  a  Trojan  orbit,  the  power  series  method 
proceeds  safely  by  tir.ie  steps  equal  to  about  1.5  canonical  units,  whereas 
a  multi-step  method  with  one  predictor  and  two  correctors  interpolating 
up  to  the  eighth  difference  had  to  progress  by  integration  intervals  4000 
times  smaller.  In  these  conditions,  the  accumulation  of  round-off  errors 
may  result  in  inaccurate  numerical  approximations  to  the  orbit  which  in 
turn  reflects  back  most  adversely  on  the  computation  of  its  fundamental 
displacements . 


II.  THE  RESTRICTED  PROBLEM 


In  the  canonical  units  of  iass,  length  and  tine,  and  in  reference 
to  the  barycentric  synodical  coordinate  system  [7],  the  planar  Restricted 
Problem  of  Three  Bodies  is  described  by  the  Hamiltonian  function 

H  m  ~  (xp  -yp  )  -  (I-iOpI1  “  UP21 

Ay  y  a 

where 

Pi  -  |  (Jrt-y)2  +  y2 
02  ■  I  (x+y-1)2  +  y2|**. 

Let  u  denote  the  4-dimensional  vector  (x,  y,  p  ,  p  ) ,  H  the  gradient 

x  y  u 

of  the  Hamiltonian  H  in  the  direction  of  the  vector  u,  H  the 

uu 

4  x  4-matrix  which  is  the  Hessian  of  H;  let  us  also  denote  by  R  a 
4  x  4-matrix  whose  elements  are  functions  of  the  time,  and  by  J  the 
4  x  4-symplectic  matrix  such  that  J  ^  ■  -J.  In  these  notations,  our 
task  can  be  described  as  the  numerical  integration  of  the  differential 
system  of  order  20 


(1) 

u  -  JH 

(2) 

R  -  JH 

the  initial  conditions  at  i  *  0  being  such  that  u(0)  does  not  lie 
either  in  the  phase  plane  (x  *  -y,  y  ■  0)  or  (x  *  1-y,  y  ■  0)  of 
binary  collisions,  and  R(0)  is  the  4  x  4-identity  matrix  I4.  The 
vector  differential  equation  (1)  yields  the  orbit,  while  the  matrix 
equation  (2)  produces  the  vesoloent  [8]  made  of  four  linearly  independent 
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variations  belonging  Co  Che  orbic. 

Because  Che  HamilConian  H  is  conservative,  Che  equaCions  admic 
Che  Jacobi  inCegral 

(3)  H  ■  consC. 
along  Che  orbic,  and  Che  vecCor  inCegral 

(4)  R*H  ■  consC. 

u 

T 

for  Che  resolvenC  of  Che  variacional  equaCions.  (N.B.  We  denote  by  R 
Che  transpose  of  the  matrix  R.) 

Because  the  equaCions  (1)  are  canonical,  the  resolvent  R  is  a 
completely  canonical  matrix,  hence  at  any  time  along  the  orbit, 

(5)  RJRT  -  J; 

this  matrix  identity  reduces  to  6  independent  bilinear  relations  between 

t 

the  fundamental  variations  along  the  orbit. 

To  the  checks  supplied  by  (3),  (4)  and  (5),  Danby  suggested  to  add 
the  vector  identity 

(6)  u(t)  -  R(t)u(0) 

which  expresses  that,  when  u  is  a  solution  of  (1),  then  u  is  a 
solution  of  the  variational  equations  (2). 

As  a  check  on  the  accuracy  of  the  numerical  integration,  the  Jacobi 
integral  (3)  is  rather  insensitive,  obviously  for  the  reason  that  it  depends 

i  » 
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not  on  the  Cartesian  components  of  the  velocity,  but  only  on  the  square 
of  the  norm  of  that  vector.  On  the  other  hand  (5)  has  proved  a  severe 
control  of  the  accuracy.  For  this  reason,  we  have  decided  to  use 

(7)  n  ■  sup  ||  R(t)JRT(t)  -  J  II* 

0<t<T 

as  a  measure  of  the  absolute  error  in  the  numerical  integration.  (N.B. 

The  norm  1 1 A |  |i  of  a  matrix  A  is  defined  as  the  sum  of  the  absolute 
values  of  its  elements.) 

We  choose  as  a  testing  stand  the  Trojan  orbit  in  the  Sun-Jupiter 

_3 

system  (M  -  0.953875  x  10  )  having  the  initial  conditions 

x  -  0.524  460  984 

y  -  0.862  013  960 

p  -  -0.850  551  063 

p  -  0.516  376  943; 

y 

this  is  an  orbit  very  close  to  being  periodic  with  period  T  •  78.505  049  481. 
Our  choice  is  dictated  by  the  fact  that  this  orbit  has  been  determined 
to  a  high  accuracy  from  another  source,  namely  from  its  representation  by 
d'Alembert  series  carried  up  to  the  fourteenth  power  of  the  orbital 
parameter  [9J. 


II.  NUMERICAL  PROCEDURES 

Among  the  available  integration  methods,  we  choose  first  the 
Runge-Kutta-Nystrom  algorithm  with  a  local  truncation  error  0(]i5), 
being  the  step-size.  A  constant  step-size  along  the  entire  period  T 
was  preferable  to  a  variable  one:  the  former  requires  only  four 


intermediate  evaluations  of  the  right  hand  members  of  (1)  and  (2) 
whereas  the  latter  imposes  six  such  evaluations  per  step  if  the  test  for 
adjusting  ii  is  to  be  performed  at  each  point.  The  error  plot  for  the 
Runge-Kutta  method  is  denoted  by  RK  in  Figure  2. 


The  second  procedure  was  the  multi-step  method  defined  by  the 
Adams-Bashf orth  predictor 


(8) 


u  - 

n-1 


+  h 


l<i<k+l 


b*z  . 
i  n-i 


and  the  Adams-Moulton  corrector 


(9) 


i  -  z  ,  +  h  7.  b.z  . 
n  n-1  i  n-i 


0<i<k 


We  based  our  decision  to  use  Adams  type  formulas  mainly  on  their  reliability. 
We  assumed  that  the  predictor-corrector  methods  for  systems  of  equations 
behave  the  same  qualitatively  as  they  do  for  single  equations;  on  this 
basis,  we  could  make  a  decision  as  to  the  type  of  methods  we  should  con¬ 
sider  [10];  our  choice  has  been  confirmed,  at  least  in  this  case  of  the 
Restricted  Problem,  by  the  experimental  results. 

We  decided  to  use  the  corrector  (9)  twice  per  step.  The  results 
showed  that  Adams-Moulton  algorithms  with  only  one  corrector  tended  to  be 
unstable,  especially  as  k  increased.  On  the  other  hand,  the  accuracy 
attained  by  three  correctors  was  only  slightly  better  than  with  two 
correctors,  for  the  same  k  and  the  same  step-size. 

The  error  behaviour  for  the  multi-step  method  with  two  corrections, 

and  k  ■  5,  6,  7,  8  successively  are  summarized  in  Figure  1.  Because 

k+2 

the  local  truncation  errors  of  formulas  (8)  and  (9)  are  0(hr1  ),  we 
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would  like  to  use  as  high  a  value  of  k  as  possible  and  still  maintain 
stability.  However ,  as  expected,  instability  appears  at  smaller  It  as 
k  increases.  Nevertheless,  since  we  are  interested  in  keeping  the  error 
measure  n  less  than  about  10  we  were  led  to  choose  eventually 
k  *  8.  The  error  plot  for  k  ■  8  is  denoted  by  PC  in  Figure  2. 


The  third  procedure  was  the  method  by  recurrent  power  series  [6]. 
Here  we  have  to  evaluate  the  coefficients  in  the  20  power  series: 


x(t+h)  -  J^X-tOh0"1, 

n^i  n 

y(t+h)  -  J^y^Oh11”1, 
n^l 

p  (t+h)  -  £p  (Oh11”1, 
X  n>l  n 


Py(t+h) 


2X 

n>l 


(t)h 


n-1 


for  the  orbit  in  its  four  dimensional  phase  space,  and 


6x(i)(t+h) 


6y(i)(t+h) 


fip^U+h) 


fip^U+h) 


Ek^coh"'1, 


E  (t)hn_1 

n 

n>l 

Euf’wh"'1, 

nil  ” 

E  v*1* (t)hn_1 

n>l  n 


for  each  column  (i  *  1,  2,  3,  4)  of  the  resolvent  matrix  R.  In  order 
to  compute  these  coefficients,  it  has  been  found  practical  to  introduce 
13  other  auxiliary  variables  to  be  computed  also  by  power  series  expansions. 
The  necessary  recurrent  formulas  can  be  found  in  [11];  they  serve  to 
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compute  all  coefficients  in  the  above  series »  starting  from  the  values 
taken  by  the  unknowns  at  time  t_Q.  Tests  are  incorporated  in  the  pro¬ 
cedure  to  decide,  with  respect  to  the  required  accuracy,  the  number  of 
terms  to  be  used  in  the  series,  and  the  step  li.  Then  the  variables 
are  computed  from  the  series  at  time  to  +  h,  and  these  new  values 
serve  to  re-initialize  the  recurrent  determination  of  the  power  series 
for  a  new  integration  step. 

In  the  problems  considered,  it  turned  out  that  the  most  practical 

procedure  was  to  keep  the  number  of  terms  in  the  series  equal  to  16 

throughout  the  entire  orbit,  and  to  vary  only  the  time  intervals  over 

-9 

which  the  series  were  employed.  To  keep  the  error  measure  n  <  10 
over  a  period,  only  about  50  points  were  necessary,  which  means  an 
unusually  large  time  step  of  1*5  unit. 

III.  COMPARISONS 

The  computations  were  performed  in  double  precision  on  an  IBM 
709^  from  programs  written  in  FORTRAN  IV. 

In  the  following  figures  we  have  plotted  for  the  various  integration 
methods  the  error  measure  n  on  the  fundamental  variations  versus  the 
computer  time  used  to  calculate  the  Trojan  orbit  over  its  period.  We 
mean  the  actual  computation  time  starting  after  the  data  were  read  in, 
and  ending  after  the  computation  was  complete,  but  before  any  printout 
occurred;  such  a  time  is  measured  iy  the  machine  clock  within  120  milli¬ 


seconds  . 
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Fig.  1 
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Fig.  2 


Such  graphs  possibly  may  be  criticized  on  the  ground  that  computer 
time  is  too  closely  dependent  on  the  programmer  and  his  ability  to  convert 
a  mathematical  procedure  into  an  efficient  sequence  of  computer  operations. 
We  point  out,  however,  that  the  Runge-Kutta  method  is  easy  enough  to 
program,  the  predictor-corrector  methods  are  of  moderately  greater 
difficulty,  whereas  the  recurrence  by  power  series  can  require  more  care. 
Thus  our  graphs,  if  anything,  are  biased  in  favor  of  the  Runge-Kutta 
and  the  multi-step  methods. 

The  general  pattern  of  the  error  plots  for  the  particular  numerical 
procedures  is  as  expected.  Basically  for  all  methods,  there  exists  a 
region — the  round-off  region — where  the  accumulation  of  rounding  errors 
is  the  main  contributing  factor  to  the  total  error;  it  lies  at  the  right 
hand  side  of  the  figures,  because  it  occurs  for  too  small  values  of  the 
step  h_,  hence  a  significant  increase  in  the  amount  of  computing  time. 
Truncation  regions  exist  at  the  intermediate  values  of  the  computing 
time;  here  the  main  contribution  to  the  total  error  is  the  accumulation  of 
the  local  truncation  errors.  In  the  round-off  region,  the  increase  in 
error  after  the  curves  passed  a  minimum  is  a  direct  result  of  the  accumu¬ 
lation  of  rounding  errors  in  the  numerical  solution  of  the  motion  equations 
(1)  rather  than  the  variational  equations  (2).  This  trend  was  indicated 
by  an  increasingly  inaccurate  Jacobi  constant  (3)  in  that  region. 

For  all  the  methods,  the  time  step  can  be  chosen  so  that  n  remains 
-9 

less  than  about  10  on  the  entire  period;  the  accuracy  can  even  be  in¬ 
creased  by  a  factor  of  100  for  the  methods  by  recurrent  power  series  and 
predictor-corrector.  In  these  conditions  of  maximum  accuracy  for  the  last 
two  methods,  the  program  returned  at  the  end  of  the  orbit  a  Jacobi  constant 


1 
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(3)  unchanged  with  respect  to  the  initial  one;  the  final  variational 
constants  (4)  were  left  constant  to  about  12  significant  figures.  The 
orbit  came  back  upon  itself  with  approximately  an  equal  number  of  signi¬ 
ficant  digits;  but  the  coordinates  at  the  end  of  the  assumed  period  T 
differed  consistently  from  the  initial  ones  by  6  figures  on  the  sixteenth 
place,  which  indicated  that  the  period  T  was  not  determined  with  sufficient 
precision. 

IV.  CONCLUSION 

From  Figure  2,  it  comes  out  quite  obviously  that  the  power  series 
is  the  most  economic  of  the  three  methods.  For  an  interval  of  n 
between  10  and  10  ^ ,  the  predictor-corrector  requires  as  much  as 
three  times  the  computing  time  employed  by  the  power  series.  In  the  same 
domain,  the  Runge-Kutta  computing  time  is  greater  than  that  of  the  power 
series  by  a  factor  between  5  and  8. 

These  conclusions  have  been  confirmed  lately  by  Dr.  Roger  Broucke 
at  the  Jet  Propulsion  Laboratory,  Pasadena,  California. 
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